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ABSTRAGT 

We use high-resolution cosmological simulations that include the effects of gasdynamics and star for¬ 
mation to investigate the origin of the Tully-Fisher relation in the standard Cold Dark Matter cosmogony. 
Stars are assumed to form in collapsing, Jeans-unstable gas clumps at a rate set by the local gas density 
and the dynamical/cooling timescale. The energetic feedback from stellar evolution is assumed to heat 
the gas surrounding regions of ongoing star formation, from where it is radiated away very rapidly. The 
star formation algorithm thus has little effect on the rate at which gas cools and collapses and, as a 
result, most galaxies form their stars very early. Luminosities are computed for each model galaxy using 
their full star formation histories and the latest spectrophotometric models. We find that at z = 0 
the stellar mass of model galaxies is proportional to the total baryonic mass within the virial radius of 
their surrounding halos. Circular velocity then correlates tightly with the total luminosity of the galaxy, 
reflecting the equivalence between mass and circular velocity of systems identified in a cosmological con¬ 
text. The slope of the relation steepens slightly from the red to the blue bandpasses, and is in fairly 
good agreement with observations. Its scatter is small, decreasing from ^ 0.45 mag in the U-band to 
^ 0.34 mag in the K-band. The particular cosmological model we explore here seems unable to account 
for the zero-point of the correlation. Model galaxies are too faint at z = 0 (by about two magnitudes) 
if the circular velocity at the edge of the luminous galaxy is used as an estimator of the rotation speed. 
The Tully-Fisher relation is brighter in the past, by about ~ 0.7 magnitudes in the B-band at z = I, 
at odds with recent observations of z ~ 1 galaxies. We conclude that the slope and tightness of the 
Tully-Fisher relation can be naturally explained in hierarchical models but that its normalization and 
evolution depend strongly on the star formation algorithm chosen and on the cosmological parameters 
that determine the universal baryon fraction and the time of assembly of galaxies of different mass. 

Subject headings: cosmology: theory - galaxies: formation, evolution - methods: numerical 


1. INTRODUCTION 

The tight correlation between the total luminosity of 
disk galaxies and the rotation speed of their gas and stars 
(Tully & Fisher 1977) has been a cornerstone of studies 
that attempt to measure galaxy distances and to map 
the dynamics of the local universe (see, e.g., Strauss & 
Willick 1995 for a review). The substantial observational 
effort devoted to the subject over the past two decades 
has provided a good understanding of the properties of 
this relation in the local universe. The Tully-Fisher (TF) 
relation steepens systematically from the blue to the red 
passbands and is surprisingly tight, especially at longer 
wavelengths (the rms scatter is only ~ 0.4 mag in the I- 
band, see Willick et al. 1997). The situation is less clear at 
higher redshift where, depending on sample selection and 
the details of the observational technique used to estimate 
rotation speeds, conflicting claims that the TF relation ei¬ 
ther brightens or dims at modest redshifts (z<l) can be 
found in the literature (Rix et al. 1997, Simard & Pritchet 
1996, Hudson et al. 1998, Vogt et al. 1996, 1997). 

From a theoretical point of view, ideas regarding the 
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origin of the TF relation can be divided into two broad 
camps: one that sees the TF relation as a result of self- 
regulated star formation in disks of different mass (see, 
e.g., Silk 1997); and another where the TF relation is a di¬ 
rect consequence of the cosmological equivalence between 
mass and circular velocity (see, e.g.. Mo, Mao & White 
1998). This equivalence is a consequence of the finite age 
of the universe, which imposes a maximum radius from 
where matter can accrete to form a galaxy. This radius 
(sometimes called the “virial” radius) corresponds roughly 
to that of a sphere with mean density about 200 times the 
critical density for closure (White & Frenk 1991, White et 
al. 1993). The circular velocity at the virial radius (V 200 ) 
and its enclosed mass (M 200 ) then are equivalent measures 
of the system mass related by 

M„. = 2.33xl0=(j^)’(^)t-Me. (1) 

H{z) is the value of Hubble’s constant at redshift z, and 
Hq = 100 h km s“^ Mpc“^ is its value at the present time. 
Given this scaling, if disk rotation speeds (Uot) are pro¬ 
portional to the circular velocity of the halo (V 200 ) and 
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stellar masses are proportional to M 200 , then luminosities 
should scale approximately like l^ot > ^ velocity dependence 
similar to that of the observed TF relation. 

Both of these assumptions may fail in practice. Disk ro¬ 
tation speeds depend in a non-trivial fashion on the stellar 
and gaseous mass of the galaxy, on the contribution of the 
dark halo within the optical radius of the galaxy, and on 
the radius where rotation speeds are measured. All these 
parameters may vary significantly from galaxy to galaxy 
(see, e.g., Navarro 1998 and references therein). Further¬ 
more, luminosities, especially in the blue passbands, are 
not simply proportional to the stellar mass, but rather re¬ 
flect a subtle interplay between the integrated and instan¬ 
taneous star formation history of each galaxy. Theoretical 
progress in this area then requires a careful evaluation of 
these effects. 

We present in this Letter a first attempt at addressing 
these issues through direct numerical simulations of the 
formation of galaxies within their full cosmological con¬ 
text. Our numerical experiments treat self-consistently 
many of the relevant physical processes of the problem: 
gravity, gas dynamics, and (in a simplified manner) the 
formation of stars and their energetic feedback. The star 
formation history of each galaxy is used as input to the 
latest spectrophotometric codes to compute broad-band 
luminosities that can be readily compared with observa¬ 
tions. Our simulations represent a significant improvement 
over prior numerical work on the subject, which either ne¬ 
glected star formation and feedback and focussed on the 
ability of gas the cool and collapse within dark halos, or 
lacked the proper numerical resolution to analyze the TF 
relation over a large range of mass and redshift (see, e.g., 
Evrard, Summers & Davis 1994, Navarro & White 1994, 
Tissera, Lambas & Abadi 1997). 

2. THE NUMERICAL EXPERIMENTS 

The simulations were performed using GRAPESPH, a 
code that combines the Smoothed Particle Hydrodynamics 
(SPH) approach to numerical hydrodynamics with a direct 
summation N-body integrator optimized for the special- 
purpose hardware GRAPE (Steinmetz 1996). GRAPE¬ 
SPH is fully Lagrangian and highly adaptive in space 
and time due to the use of individual particle smooth¬ 
ing lengths and timesteps. It is thus optimally suited to 
study the formation of highly non-linear systems such as 
individual galaxy systems in a cosmological context. The 
code used for the simulations described in this paper in¬ 
clude the self-gravity of gas, stars, and dark matter, a full 
3D hydrodynamical treatment of the gas, radiative and 
Gompton cooling, and a simple recipe for transforming 
gas into stars. 

The star formation algorithm we adopt is similar to that 
described in Steinmetz & Muller (1994, 1995, see also Katz 
1992 and Navarro & White 1993). Star formation is mod¬ 
eled by creating new collisionless “star” particles in col¬ 
lapsing regions that are locally Jeans-unstable at a rate 
given by p* = c*pgas/max(Tcooi, Tdyn)- Here Pgas is the 
gas density and TcooI and Tdyn are the local cooling and 
dynamical timescales, respectively. The proportionality 
parameter, c* = 0.05, is chosen so that most eligible gas 
is transformed into stars in about ~ 20 Tdyn in dense re¬ 
gions, where TcooI ^ Riyn- Our results at z = 0 are weakly 
dependent on our choice of c* since dynamical and cool¬ 


ing timescales in regions that actively form stars are much 
shorter than the age of the universe. 

After formation, “star” particles are only affected by 
gravitational forces, but they devolve energy to their sur¬ 
roundings in a crude attempt to mimic the energetic feed¬ 
back from evolving stars and supernovae. Star particles 
inject 10^® ergs per solar mass of stars formed into their 
surrounding gas 10^ yrs after their formation. This energy 
is invested into raising the internal energy (temperature) 
of the gas. However, because stars form in high-density 
regions, where cooling timescales are short, the feedback 
energy is almost immediately radiated away, and as a re¬ 
sult the star formation history traces the rate at which 
gas cools and collapses to the center of dark matter halos. 
Our implementation may thus be thought of as a “minimal 
feedback” model. We caution that other possible imple¬ 
mentations, in particular those that include modifications 
to the kinetic energy of the gas surrounding star forming 
regions, may lead to quite different star formation histories 
than the ones we report here (Navarro & White 1993). 

The cosmological model we investigate is the standard 
Gold Dark Matter (CDM) scenario, the paradigm of hier¬ 
archical clustering models of structure formation, with the 
following choice for the cosmological parameters: D = 1, 
h = 0.5, Hb = 0.0125 and A = 0. The power spec¬ 
trum is normalized so that at z = 0 the rms amplitude 
of mass fluctuations in 8h“^Mpc spheres is as = 0.63. 
Although long discredited as a viable model for structure 
formation in the universe, the standard Gold Dark Matter 
scenario is the best studied hierarchical clustering model 
and therefore it serves as a useful cosmological testbed of 
different galaxy formation models. 

We simulate spherical regions that evolve to form ha¬ 
los with circular velocities in the range (80, 350) km s“^. 
These regions are identified in cosmological simulations of 
large periodic boxes and resimulated individually at higher 
resolution, including the tidal field of the original simula¬ 
tion. Each region is filled initially with ~ 32, 000 gas par¬ 
ticles and the same number of dark matter particles. The 
size of the sphere is scaled in each case so that halos have 
about the same number of particles at z = 0. This en¬ 
sures comparable numerical resolution in systems of very 
different mass. Details of the procedure can be found in 
Navarro & White (1994). The gas particle mass varies 
between 2.5 x 10® and 5 x IO^/i^^Mq, depending on the 
mass of the system under consideration. Dark matter par¬ 
ticle masses scale in the same way and vary from 5 x 10^ 
to 10^h~^M q. These particle masses are small enough to 
exclude artificial heating effects due to collisions between 
gas and dark matter particles (Steinmetz & White 1997). 
All simulations start at z = 21 and use gravitational soft¬ 
enings between 0.5 and 1.25h~^ kpc. 

Stars form preferentially in collapsed, dense clumps of 
gas, and therefore model “galaxies” are easily identified 
in these simulations as dense collections of star parti¬ 
cles. In the interest of simplicity, we retain for analysis 
only the most massive galaxy inhabiting each individual 
halo. Most stars in a galaxy are contained within a radius 
i"gai = 15(V2oo/ 220 km s~^) h~^ kpc, which we use as a 
fiducial radius to compute all properties of the luminous 
component of model galaxies. We also cull from the sam¬ 
ple all ongoing major mergers, defined as galaxies where 
the rms velocity of the stars differs from the circular ve- 
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locity at rgai by more than 30%. In practice less than 10% 
of the systems are eliminated this way from the analysis. 
The final sample consists of 93 galaxies at z = 0 and of 44 
galaxies at z = 1, most of which contain more than 1,000 
star particles. 

Galaxy luminosities are computed by simply adding up 
the luminosities of each star particle, taking into account 
the time of creation of each particle (its “age”) and its 
metallicity, as described in detail by Contardo, Steinmetz 
& Fritze-von Alvensleben (1998). The latest version of the 
spectrophotometric models of Bruzual & Chariot (1993) is 
also used for comparison. 

3. RESULTS AND DISCUSSION 
3.1. The Tband Tully-Fisher relation at z = 0 

Figure 1 shows the simulated Tully-Fisher relation at 
z = 0 in the I-band, compared with the data of Gio- 
vanelli et al. (1997), Mathewson, Ford & Buchhorn (1992) 
and Han & Mould (1992). The lower panel shows I- 
band magnitudes versus V 200 , the circular velocity at the 
virial radius; the upper panel uses instead Vrot, the cir¬ 
cular velocity measured at the fiducial radius rgai. The 
slope and scatter of the simulated TF relation are in fairly 
good agreement with the observational data, regardless of 
which velocity estimator is used. This result also holds in 
other bandpasses: the model TF relation becomes shal¬ 
lower (and the scatter increases) towards the blue, just as 
in observational samples (see Table 1). 

The zero-point is also in good agreement with the data, 
provided that V 200 is used to estimate rotation speeds. 
Circular velocities at rgai (the edge of the luminous galax¬ 
ies) are about 40-60% higher, so that adopting Vrot as a 
velocity estimator results in a TF relation almost two mag¬ 
nitudes too faint at given rotation speed. The large differ¬ 
ence between Vj-ot and V 200 is caused by the large fraction 
of baryons that assemble into the central galaxy and the 
dark mass they draw in as they settle at the center of the 
halo. On average, about 70-90% of all baryons within the 
virial radius are confined within rgai at z = 0. 


Band 

Mo - 5 \og{h) 

5 Mo 

A 


cr(M) 

U 

-18.50 

0.066 

6.95 

0.30 

0.45 

B 

-18.66 

0.060 

7.08 

0.28 

0.39 

V 

-19.34 

0.058 

7.20 

0.25 

0.35 

R 

-19.89 

0.056 

7.30 

0.25 

0.34 

I 

-20.27 

0.060 

7.33 

0.25 

0.34 

K 

-21.95 

0.053 

7.48 

0.23 

0.34 

Table 1: 

Parameters of 

fits of the form M 

= Mo + 


Alog(Vrot/200kms“^) to the TF relation in different bands. 
Errors in Mq and A are (10,90) percentile bootstrap devia¬ 
tions from the mean. a{M) is the rms scatter in magnitudes 
of a least square fit. 

This illustrates what is likely to be a fatal problem for 
this particular cosmological model. Because dark halos 
formed in this scenario are quite centrally concentrated 
(Navarro et al. 1997), the assembly of a massive galaxy at 
the center raises Vj-ot above and beyond the halo circular 
velocity, by up to 60% in the standard CDM scenario we 
investigate here. The effect of this increase is a significant 
mismatch of the zero-point in the model TF relation rel¬ 
ative to observations. Model galaxies would have to be 


almost 2 magnitudes brighter in order to match the ob¬ 
served TF relation. This is highly unlikely, since galaxies 
in our model contain most of the baryonic material within 
the virial radius. 

The only way to collect a massive disk galaxy without 
significantly increasing Vrot over V 200 is to have dark ha¬ 
los which are less dense than those formed in the standard 
CDM scenario. This favors low-density (Dq < 1) universes, 
where dark halos are less centrally concentrated (Navarro 
et al 1997, Navarro 1998). In this case, the problem is 
further alleviated because the universal baryon fraction 
is larger (by a factor Dq ^) and therefore massive, bright 
galaxies can be assembled by collecting a relatively small 
fraction of baryons into the central galaxy. 
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Fig. 1.—I-band Tully Fisher relation at z = 0. The error 
bars in the simulated data indicate the difference in magni¬ 
tudes that results from adopting a Salpeter or a Scalo IMF. 


3.2. The star formation history of model galaxies 

As discussed in §2, our feedback implementation has 
minimal impact on the rate at which gas cools and col¬ 
lapses to the center of dark matter halos. The star form¬ 
ing history of a model galaxy thus traces the history of 
cooling and accretion of gas within halos, punctuated by 
small “bursts” of star formation that coincide with ma¬ 
jor mergers. This is shown in Figure 2, where we plot the 
combined star formation rate of four model galaxies, span¬ 
ning a factor of ^ 3 in circular velocity. The large peak 
at early times is due to the nearly simultaneous collapse 
of a number of halos which later merge to form the final 
galaxy. The star formation rate in each progenitor rarely 
exceeds ~ 5OM0/yr, although its combined rate can be 
much higher. 

Most simulated galaxies form their stars quite early (half 
of them typically by z ^ 1-1-5), and have almost exhausted 
their gas supply at z = 0, when the star formation rate 
is only about 20-30% of the global rate averaged over the 
age of the universe. These star formation histories are thus 
similar to those deduced for E/SO/Sa galaxies rather than 
to the nearly constant star formation rates characteristic 
of late-type spirals (Kennicutt 1998). 
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Fig. 2.—Star formation rate smoothed over 10® yr intervals 
(in solar masses per year) versus lookback time (in Gyr) for 4 
typical halos with virial velocities between 90 and 250 km s“ . 


3.3. The angular momentum of simulated galaxies 

As discussed in the previous subsection, most gas cools 
and collapses before the merger events through which the 
galaxy is finally assembled. We therefore expect the an¬ 
gular momentum of the model galaxies to be significantly 
lower than those of disk galaxies, since angular momentum 
is transferred from the gas to the dark matter during merg¬ 
ers (Navarro, Frenk & White 1995, Navarro & Steinmetz 
1997). 

This is shown in Figure 3, where we compare the spe¬ 
cific angular momentum of the model galaxies to the disk 
galaxies in the samples of Courteau (1997), Mathewson et 
al. (1992), and the compilation of Navarro (1998). Disk 
angular momenta are estimated assuming that they are 
pure exponential systems with flat rotation curves, jdisk = 
2 Vrot ?’disk, where K-ot is identified with the (maximum) ro¬ 
tation speed and ruisk is the exponential scalelength. The 
solid line shows, for comparison, the typical specific angu¬ 
lar momentum of dark halos, jhaio = (2//c)^^^ A (V 2 oo/km 
s“^)^ km s~^h~^ kpc, where fc is a factor of order unity 
that depends on the internal structure of halos (Mo et 
al. 1997) and A is the usual dimensionless spin parame¬ 
ter, typically found to be ~ 0.05 in cosmological N-body 
simulations (see, e.g.. Cole & Lacey 1996). (We assume 
k 2 oo=Kot to plot the halo curve.) Disk angular momenta 
scale with rotation speed just as jhaio scales with circular 
velocity, suggesting that baryons in disks have retained 
a roughly constant fraction of the angular momentum of 
their surrounding halos: about one half if V 200 = Vrot (see 
Syer, Mo & Mao 1998 for a more detailed discussion). 

The simulated galaxies, on the other hand, have spe¬ 
cific angular momenta which are more than one order of 
magnitude lower than those of observed disks. This con¬ 
firms that the hierarchical assembly of a galaxy robs the 
baryonic component of much of its angular momentum if 
gas cooling and star formation precedes most major merg¬ 
ers. This angular momentum problem signal that feedback 
should be more effective at preventing gas from collapsing 
and cooling in the early stages of the hierarchy in order to 
ensure that most of the angular momentum is not lost as 
the luminous component of a galaxy is assembled (see also 
Weil, Eke & Efstathiou 1998). 



Fig. 3.—Specific angular momentum vs circular velocity of 
model galaxies compared to observational data. 

3.4. The redshift evolution of the Tully-Fisher relation 


Figure 4 shows the redshift evolution of the TF relation 
in the rest-frame B-band. The shaded region corresponds 
to the results at z = 0, as reported by Pierce & Tully 
(1992). Starred symbols show the z ~ 1 data of Vogt et 
al. (1996, 1997). All data have been scaled to the value 
of go = 0-5 used in this paper. Note that for this value of 
go the data from Vogt et al. are consistent with a slight 
dimming of the TF relation at z ~ 1. Their claim of a mod¬ 
est brightening holds only for the low value of go = 0.05 
adopted in their analysis. 

At z = 0 the results for the B-band TF relation of simu¬ 
lated galaxies are similar to those presented in Figure 1 for 
the I-band. The slope and scatter are in good agreement 
with observations, but the zero-point depends strongly on 
the radius at which circular velocities are measured. Good 
agreement is found provided that V 200 is used as a measure 
of the rotation speed, but a significant zero-point offset is 
observed if Vot is used instead. 

Independent of which velocity estimator is adopted, 
the model TF relation brightens at z = 1, by ~ 0.7 or 
~ 0.2 mag, depending on whether Vot or V 200 is used, 
respectively. This is at odds with the observational re¬ 
sults and perhaps a bit surprising since masses scale like 
oc (1 -I- at fixed circular velocity (eq. 1) and 

are, therefore, a factor of ~ 2.8 times smaller at z = 1. 
However, star formation rates are much higher at z = 1 
(by a factor of ~ 4) and the decrease in mass is more 
than compensated by the higher abundance of young stars 
present at that time. 

This interpretation is confirmed by results in other pass- 
bands less sensitive to the instantaneous star formation 
rate: the K-band TF relation actually dims by about 0.1- 
0.4 mag between z = 0 and 1 (depending again on whether 
V 200 or Vrot is used). The brightening trend observed in 
Figure 4 is therefore a direct consequence of the star for¬ 
mation algorithm we use, and not a genuine feature of 
hierarchical clustering models. One example is provided 
by the study of Mao, Mo & White (1998), who use dif¬ 
ferent assumptions for the star formation history of disks, 
and argue for a roughly constant TF zero-point between 
z = 0 and 1 in a cosmological model similar to the one we 
adopt here. 
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Fig. 4.—Redshift evolution of the rest-frame B-band Tully- 
Fisher relation. 


4. SUMMARY 

We have used high-resolution cosmological simulations 
which incorporate the effects of gas pressure, hydrodynam- 
ical shocks, radiative cooling, star formation and super¬ 
nova feedback to study the origin of the TF relation of 
spiral galaxies. Our main conclusion is that the slope (and 
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scatter) of this relation is naturally reproduced in hierar¬ 
chically clustering universes, and reflect the equivalence 
between mass and circular velocity of systems formed in 
a cosmological context. The particular model we explore 
here (the standard CDM scenario) seems unable to ac¬ 
count for the zero-point of the correlation. Model galax¬ 
ies are too faint at z = 0 (by about two magnitudes) if 
the circular velocity at the edge of the luminous galaxy 
is used as an estimator of the rotation speed. This re¬ 
sult had already been hinted at, although using different 
arguments, by semianalytic models of galaxy formation 
(Kauffmann, White & Guiderdoni 1993, Cole et al. 1994, 
but see Somerville & Primack 1997 for a dissenting view). 

The model TF relation brightens in the past, at odds 
with current observational evidence. However, this is a 
direct result of the particular star formation algorithm we 
use, which transforms gas into stars very efficiently at high 
redshift. Efficient gas cooling (and star formation) that 


precedes the final assembly of the galaxy through mergers 
leads to the formation of stellar systems with angular mo¬ 
menta far lower than those of observed disks. The combi¬ 
nation of successes and failures we report here demonstrate 
the potential of galaxy scaling laws and their evolution as 
probes of the intricate physics of star formation and feed¬ 
back operating in a cosmological context. 
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